Spin Stiffness of Stacked Triangular Antiferromagnets 
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We study the spin stiffness of stacked triangular antiferromagnets using both heat bath and broad 
histogram Monte Carlo methods. Our results are consistent with a continuous transition belonging 
to the chiral universality class first proposed by Kawamura. 
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I. INTRODUCTION 

The magnetic ordering of geometrically frustrated an- 
tiferromagnets differs substantially from the convai 
tional magnetic ordering in nonfrustrated magnetstH 
Indeed, the nature of the phase transition in the 
case of starljprj n'T'if .TE-'.'^y antiferromagnets remains 
controversiaffleml§t3BEl. The geometry of the 
stacked triangular lattice has triangles as the elemen- 
tary units and this arrangement inhibits an anti-parallel 
alignment of the spins in each triangular layer. Conse- 
quently, the system is said to be geometrically frustrated. 
This frustration leads to a comprise where the spins on 
each triangle adopt a non-coUinear spin ordering at low 
temperatures. The spins form a planar configuration in 
which nearest neighbours are oriented at an angle of 120 
degrees with respect to one another. The ground state 
can be described by a matrix like order parameter giving 
the orientation of each spin on the elementary triangles 
and forms an SO (3) parameter spacetJ. This unusual 
symmetry of the order parameter and the appearance 
of 'chiral' degrees of freedom which correspond to the 
ground state having two possible realizations, left and 
right handed, lead KawamuraEJ to conjecture the exis- 
tence of a new chiral universality class . The chiral de- 
grees of freedom are believed to be responsible for the 
novel critical behaviour but they are not decoupled from 
the spin degrees of freedom and the two quantities order 
simultaneously. While recent field-theoretic renormaliza- 
tion group studies of this system using an expansion up 
to six loops in fixed dimension d — i indicate the ex- 
istence of a stable fixed point that |-aorresponds to the 



proposed chiral universality clasalU'O, similar analyses 
using a three loop perturbation technique as well as an 
epsilon expansion approach to the same order, did not 
find a stable fixed point and hence exclude the possibil- 
ity of ai-continuous phase transition for this frustrated 
systemQEi. Non perturbative RG approaches find that 
the phase transition is possibly a very weak first order 
transition with effective critical exponentaJ. 

In the present work we use both a standard heat 
bath Monte Carlo method as skeU as a recently devel- 
oped broad histogram methodEj to study the classical 
isotropic antiferromagnet on this geometry. In particu- 



lar, we study the spin stiffness which provides a direct 
measure of the correlation length exponent v. The spin 
stiffness is a convenient quantity to study since it does 
not require knowledge of the order parameter but does 
measure the rigidity of the ordered phase against fluctu- 
ations. Our results confirm the picture of a continuous 
transition belonging to a new chiral universality class. 



II. MODEL AND METHODS 

The model is described by the following Hamiltonian 
7J=-^JS',-5, -^J'5fc-5,. (1) 
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where Si is a classical three component vector of unit 
length located at the sites i of a hexagonal lattice. The 
first sum is over nearest neighbours in the triangular 
planes which interact with an antiferromagnetic coupling 
J < and the second sum is over inter-plane nearest 
neighbours which are taken to have a ferromagnetic cou- 
pling J' > with I J' I = I J| = 1. Hence all energies and 
temperatures are measured in units of | J|. 

We study the response of the system to a virtual 
twist of |the spin system. The spin stiffness , or helicity 
modulusli^ , measures the increase in free energy associ- 
ated with twisting the order parameter in spin space by 
imposing a gradient of the twist angle about some axis n 
in spin space along some direction u in the lattice. The 
spin stiffness can be written as a second derivative of 
the free energy with respect to the strength of the gra- 
dient arid, can be calculated as an equilibrium response 
functiontZl . Finite size scaling theory predicts that the 
spin stiffness should vanish at the critical point with an 
exponent related to the correlation length exponent. 

We calculate the diagonal elements of the spin stiffness 
tensor corresponding to twists about three principal di- 
rections in spin space. If we divide the lattice sites into 
three equivalent sublattices A, B and C corresponding to 
the vertices of the elementary triangles, then a chirality 
vector can be defined to characterize the non-coUinear 
ordering of the spins. The chirality is defined locally for 



2 



each upward (downward) triangle by the following ex- 
pression 

Ka^SaxSe+SbX-Sc + Scx Sa (2) 

In the ground state, the chirality is uniform and per- 
pendicular to the spin planes. This symmetry of the 
order parameter suggests that the average chirality di- 



rection {K) be chosen for one of the principal axes and 
the other two directions (li, J_2) are chosen to be in the 
spin plane perpendicular to the average chirality vector 
such that the three vectors form an orthonormal triad. 
The spin stifiaess component Pa at temperature T can 
be written asE-3 



1<J 
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(3) 



where a,/3,7 = if, ±i,_L2 and the indices are taken in 
cyclic order. The twist is taken to be along the u direc- 
tion in the lattice and e-ij is a unit vector directed along 
the nearest neighbour bond from site i to j . The angu- 
lar brackets indicate a thermal average in the canonical 
ensemble. Since the ground state is a planar spin arrange- 
ment, the stiffnesses satisfy a perpendicular axis theorem 
Pj^ = Pjl^ + Pj_^ at zero temperature. Deviations from 
this relationship are a measure of fluctuations of spins 
from the planar order. 

We perform numerical simulations using both a con- 
ventional Monte Carlo (MC) heat bath method and the 
more recent broad histogram method (BHM) introduced 
by Oliveira et. al. The latter method is based on the mi- 
crocanonical ensemble approach to statistical sampling 
at fixed energy and allo wa an . accurate estimate of the 
energy density of statest3ll§ll3 g{E) . By knowing the 
density of states g{E) and the microcanonical averages 
of various quantities < Q >e, their temperature depen- 
dence can be determined by using the following expres- 
sion for the canonical averages 
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FIG. 1: Natural logarithm of the energy density of states 
versus the energy per site e for a 42 x 42 x 42 lattice in the 
range —2.1 < e < —0.5. 



III. RESULTS 
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Ee<Q >E9{E)e 



-E/T 



(4) 



In the conventional heat bath method temperature is 
tuned as an external parameter and number of tempera- 
ture points is limited by number of computer runs. The 
BHM method allows us to probe the system in a continu- 
ous range of T but requires a large number of energy bins 
for large system sizes. We simulate spin systems of size 
iV = ^iti^ ^ ^ 24, 30, 42, 60 and 66 for the heat bath 
method and only up to L = 60 for the BHM method. Pe- 
riodic boundary conditions are employed for both meth- 
ods. We find excellent agreement between these two nu- 
merical methods. 



The broad histogram method(BHM) is based on the 
relation 

g{E) < N,p{E) >= g{E + AE) < Ndn{E + AE) > (5) 

where < Nup{E) >, < Ndn{E) > are microcanonical av- 
erages which measure the number of moves which in- 
crease (decrease) the energy by the amount AE'. Once 
these microcanonical averages are known, the micro- 
canonical temperature Tm{E) can be determined from 



l/T,r,{E) = 



d\ng{E) 
dE 



:ln 



< Nup{E) > 



AE <NdniE + AE)> 



(6) 



and we can then integrate this expression in some range of 
energy to obtain the energy density of states In g{E) as a 
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FIG. 2: Energy per site e and specific heat Cv obtained using 
the BHM method for L = 24, 30, 42, 60. 



function of E. In our case the energy is a continuous vari- 
able and we divide the energy axis into bins of a fixed size 
IS.E = 1.8 such that Ai? « E, where E is total energy 
of interest. We employ a simple microcanonical dynamics 
to sample phase space and the energy density of states 
g{E) (up to a multiplicative constant) is determined us- 
ing the BHM relation above. One microcanonical sweep 
consists of a random sweep through the lattice sites and 
generating a new configuration of the spins by restricting 
the choice of a new random orientation of the spin at site 
i with respect to the local field of the nearest neighbours 
such that the total energy of the system remains within 
the energy interval E,E + AE. At any given value of E 
, 75 microcanonical sweeps were performed and 25 sam- 
ple measurements were taken of various thermodynamic 
quantities such as the energy, specific heat and spin stiff- 
ness. Before sampling the next energy interval, 40 initial 
microcanonical sweeps were performed to avoid correla- 
tions. This procedure was repeated using different seeds 
for random numbers and errors were determined using 
the standard deviations for these separate measurements. 
Figure 1 shows our results for hig{e) as a function of 



FIG. 3: Spin stiffnesses as a function of T. a) the points 
indicate the heat bath results and the lines correspond to the 
BHM results. All three stiffnesses vanish at the same finite 
temperature near T ~ .96. b) the heat bath results for pK in 
a smaller temperature range show significant finite size effects 
near Tc- 



e = E /N in the case ofa42x42x42 lattice. The units are 
arbitrary since we integrate equation (6) starting from 
e = —2.1 and not the ground state value Cq = —2.5 .The 
number of energy bins used for this energy range was 
61740. For general values of L, the number of energy bins 
required to study this same range with the same fixed size 
of energy bin is /6 and is thus of the same order as 
the number of sites. When the energy density of states 
is combined with the microcanonical averages < Q >e 
for various thermodynamic quantities, we can then plot 
them as continuous functions of T using equation (4). 
Figure 2 shows the energy per site and the specific heat 
obtained using the BHM method for various linear sizes 
L. The energy displays strong finite size effects near the 
temperature where the specific heat has a maximum. The 
figures clearly indicate that a phase transition occujb near 
T 0.96 in agreement with previous MC studies.El 

We have used both the BHM method and a Monte 
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FIG. 4: The same data as in figure 3(b) is plotted as a function FIG. 5; Lpx versus T for various lattice sizes are indicated by 

of 1/1/ for a set of equally spaced temperatures in the range the points. The lines are linear interpolations which indicate 

.85 < T < .95. Extrapolation to the large L limit yields a unique crossing point at Tc = 0.958 ± 0.002 
estimates for pK for an infinite lattice. 



Carlo heat bath method at fixed values of T to calcu- 
late the spin stiffness. In the heat bath method, we dis- 
card the first 1000 sweeps and perform 45000 MC steps 
in each run. Figure 3(a) shows both our heat bath re- 
sults, indicated by points, and the BHM results, indi- 
cated by lines, for the three stiffnesses for various lattice 
sizes L as a function of the temperature T. The relation 
Pk ~ Pli + PL2 ^^^^ satisfied for all values of T < 0.95 
indicating that there is a relatively small deviation from 
the planar spin configuration. All three stiffnesses are 
nonzero at low T and vanish near T .96 which corre- 
sponds to the specific heat divergence in figure 2. Figure 
3(b) shows the heat bath data for on an enlarged 
temperature scale. The stiffnesses clearly show large fi- 
nite size effects and approach zero near T ^ .96. The 
points labelled infinity are obtained by plotting ver- 
sus 1/L at various values of T and extrapolating to the 
large L limit as shown in figure 4. 

These finite size effects can be used to determine the 
correlation length exponent v directly. Finite size scaling 
considerations for p{T, L) predict 
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p(r,L) = -/(i/0--/(i^/ni|) 



(7) 



where t is the reduced temperature. This form suggests 
that we can plot Lp(T, L) versus T to identify Tc as the 
temperature where the curves for different values of L in- 
tersect. Figure 5 shows our heat bath results for Lp^ as a 
function of T for lattice sizes L — 24, 30, 42, 60, 66. Linear 
interpolations of ncigbouring temperature points indicate 
that the curves intersect at a value of Tc — 0.958 ±0.002. 
We have also used our BHM results in the same temper- 
ature range and we obtain the same estimate for Tc- 

In the limit as L 00, the scaling form predicts 
p ~ jij'^. Using the values of the stiffness obtained by 
extrapolating to large values of L as in figure 4 and then 
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FIG. 6: A In-ln plot of pK as L 00 versus \t\ using the 
estimated value of Tc yields a value oi v = .589 ± .007. 



plotting these versus \t\ on a In-ln scale, we can obtain an 
estimate of v. Figure 6 shows our results for pK which 
yields the value v = .589 ± .007. This value agrees very 
well with previous Monte Carlo estimatesEJ but is slightly 
larger than the value found by the recent six loop .re.pqr- 
malization group calculations in three dimensions .E3Ej 

Figure 7 shows a finite size scaling plot of our stiffness 
results using the values of Tc and v quoted above. The 
data obtained from both the heat bath MC method for 
sizes L = 24, 30, 42, 60, 66 and the BHM method for L = 
24,30,42 collapse very well to a universal function for 
temperatures below Tc- The value v — .589 ± .007 is 
certainly very different from the value v = 0.7113 which 
describjes the three dimensional Heisenberg universality 
class.Ea 
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IV. SUMMARY 

We have calculated the spin stiffness of the isotropic 
Heisenberg antiferromagnet on the stacked triangular ge- 
ometry using both a MC heat bath and BHM method. 
The spin stiffness has the advantage that it measures the 
rigidity of the ordered phase in response to a virtual twist 
without having to specify the order parameter. The re- 
sults obtained from both numerical approaches agree and 
predict a continuous phase transition which belongs to 
the new chiral universality class proposed by Kawamura. 
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